#include "Solution.h"
using namespace blitz;

bool Solution::Analysis()
{
	Initial();
	ofstream output("current.dat");
	for (int j=1; j<1000; j++)
	{
		for (int m=0; m<Object->GetEdgeNum(); m++)
		{
			emx::Real3D E = Ex->GetIncident(Object->GetEdge(m)->MidPoint,j-1);
			emx::Real3D rho = Object->GetEdge(m)->Rho_plus + Object->GetEdge(m)->Rho_minus;
			emx::real_data rhs = dot(rho,E)*Delta - dot(rho,_A(m,j)-A(m,j-1)) + 2 * Delta * (Phi(Object->GetEdge(m)->LeftFacet->Idx,j-1)-Phi(Object->GetEdge(m)->RightFacet->Idx,j-1));
			I(m,j) = rhs / dot(K(m,m),rho);
			if (j > 2) I(m,j-1) = (I(m,j-2)+2*I(m,j-1)+I(m,j))/4.0;
		}
		output << I(42,j) << std::endl;
	}
	output.close();
// 	I(0,1) = 9.125;
// 	I(0,2) = 2.23123;
// 	std::cout << I_Integral(0,1.5*Delta) << std::endl;
// 	std::cout << 0.875*Delta*I(0,1) + 0.125*Delta*I(0,2) << std::endl;
	return true;
}

emx::real_data Solution::T(int j, emx::real_data t)
{
	if ((j-1)*Delta>=t || (j+1)*Delta<=t) return 0;
	else return 1-fabs(t-j*Delta)/Delta;
}

bool Solution::Initial()
{
	K.resize(Object->GetEdgeNum(),Object->GetEdgeNum());
	phi.resize(Object->GetPatchNum(),Object->GetEdgeNum(),2);
	I.resize(Object->GetEdgeNum(),1210);
	for (int m=0; m<Object->GetEdgeNum(); m++){
		for (int k=0; k<Object->GetEdgeNum(); k++){
			K(m,k) = Calculate_K(m,k);
		}
		I = 0.0;
	}
	for (int m=0; m<Object->GetPatchNum(); m++){
		for (int k=0; k<Object->GetEdgeNum(); k++){
			phi(m,k,0) = Calculate_phi(m,k,0);
			phi(m,k,1) = Calculate_phi(m,k,1);
		}
	}
	Delta = Ex->GetDelta();
	return true;
}

emx::real_data Solution::Singular(int patchIdx)
{
	emx::Real3D r[3], rho[3];
	emx::real_data p[3], sinBeta[3], _sinBeta[3], cosBeta[3], _cosBeta[3];
	emx::real_data result = 0.0;
	for (int i=0; i<3; i++)
	{
		r[i] = Object->GetPatch(patchIdx)->Vertex((i+1)%3)->Position - Object->GetPatch(patchIdx)->Vertex(i)->Position;
		rho[i] = Object->GetPatch(patchIdx)->Vertex(i)->Position - Object->GetPatch(patchIdx)->Center;
		sinBeta[i] = -dot(rho[i],r[i])/(sqrt(dot(rho[i],rho[i])*dot(r[i],r[i])));
		cosBeta[i] = sqrt(1-sinBeta[i]*sinBeta[i]);
		p[i] = cosBeta[i]*sqrt(dot(rho[i],rho[i]));
		_sinBeta[i] = dot(rho[(i+1)%3],r[i])/(sqrt(dot(rho[(i+1)%3],rho[(i+1)%3])*dot(r[i],r[i])));
		_cosBeta[i] = sqrt(1-_sinBeta[i]*_sinBeta[i]);
	}
	for (int i=0; i<3; i++)
	{
		result += p[i]*log(cosBeta[i]*(1+_sinBeta[i])/(_cosBeta[i]*(1-sinBeta[i])));
	}

	return result;
}

emx::Real3D Solution::Calculate_K(int m, int k)
{
	emx::Real3D mid = Object->GetEdge(m)->MidPoint;
	emx::Real3D value(0,0,0);
	for (int i=0; i<4; i++)
	{
		emx::real_data R_plus = sqrt(dot(mid-Object->GetEdge(k)->LeftFacet->Point(i),mid-Object->GetEdge(k)->LeftFacet->Point(i)));
		emx::Real3D Rho_plus = Object->GetEdge(k)->LeftFacet->Point(i) - Object->GetEdge(k)->LeftPoint;
		emx::real_data R_minus = sqrt(dot(mid-Object->GetEdge(k)->RightFacet->Point(i),mid-Object->GetEdge(k)->RightFacet->Point(i)));
		emx::Real3D Rho_minus = Object->GetEdge(k)->RightPoint - Object->GetEdge(k)->RightFacet->Point(i);
		value += 0.25 * (Rho_plus/R_plus + Rho_minus/R_minus);
	}
	return emx::Mu0/(4*emx::PI)*(Object->GetEdge(k)->EdgeLength/2.0) * value;
}

emx::real_data Solution::Calculate_phi(int m, int k, int flag)
{
	int sourceIdx, fieldIdx = m;
	emx::real_data sign, result = 0.0;
	if (0 == flag){
		sourceIdx = Object->GetEdge(k)->LeftFacet->Idx;
		sign = -1.0;
	}
	else {
		sourceIdx = Object->GetEdge(k)->RightFacet->Idx;
		sign = 1.0;
	}

	if (sourceIdx != fieldIdx) 
	{
		for (int j=0; j<4; j++)
		{
			result += 0.25/sqrt(dot(Object->GetPatch(sourceIdx)->Point(j)-Object->GetPatch(fieldIdx)->Center,Object->GetPatch(sourceIdx)->Point(j)-Object->GetPatch(fieldIdx)->Center));
		}
	}
	else
	{
		result += 1.0/Object->GetPatch(sourceIdx)->Area*Singular(sourceIdx);
	}

	return sign*result*(Object->GetEdge(k)->EdgeLength/(4*emx::PI*emx::Epsilon0));
}

emx::real_data Solution::I_Integral(int edgeIdx, emx::real_data t)
{
	if (t<0) return 0;
	int q = int(t/Delta) - 1;
	emx::real_data value = Delta * I(edgeIdx,q+1) / 2.0;
	for (int j=1; j<=q; j++)
	{
		value += I(edgeIdx,j) * Delta;
	}
	value += (1- ((q+2)-t/Delta)*((q+2)-t/Delta))*(Delta * I(edgeIdx,q+1) / 2.0);
	value += (t/Delta-(q+1))*(t/Delta-(q+1))*(Delta * I(edgeIdx,q+2) / 2.0);
	
	return value;
}

emx::real_data Solution::Phi(int patchIdx, int j)
{
	emx::real_data value = 0.0;
	for (int k=0; k<Object->GetEdgeNum(); k++)
	{
		emx::real_data Rk_plus = sqrt(dot(Object->GetEdge(k)->LeftFacet->Center - Object->GetPatch(patchIdx)->Center,Object->GetEdge(k)->LeftFacet->Center - Object->GetPatch(patchIdx)->Center));
		emx::real_data tR_plus = j*Delta - Rk_plus/emx::C;
		emx::real_data Rk_minus = sqrt(dot(Object->GetEdge(k)->RightFacet->Center - Object->GetPatch(patchIdx)->Center,Object->GetEdge(k)->RightFacet->Center - Object->GetPatch(patchIdx)->Center));
		emx::real_data tR_minus = j*Delta - Rk_minus/emx::C;
		value += I_Integral(k,tR_plus) * phi(patchIdx,k,0);
		value += I_Integral(k,tR_minus) * phi(patchIdx,k,1);
	}
	
	return value;
}

emx::Real3D Solution::A(int edgeIdx, int j)
{
	emx::Real3D value(0,0,0);
	for (int k=0; k<Object->GetEdgeNum(); k++)
	{
		emx::real_data Rmk = sqrt(dot(Object->GetEdge(edgeIdx)->MidPoint - Object->GetEdge(k)->MidPoint,Object->GetEdge(edgeIdx)->MidPoint - Object->GetEdge(k)->MidPoint));
		emx::real_data t = j*Delta - Rmk/emx::C;
		if (t<0) continue;
		int floor = int(t/Delta);
		int ceiling = floor + 1;
		value += (T(floor, t) * I(k,floor) + T(ceiling, t) * I(k,ceiling)) * K(edgeIdx,k);
	}

	return value;
}

emx::Real3D Solution::_A(int edgeIdx, int j)
{
	emx::Real3D value(0,0,0);
	for (int k=0; k<Object->GetEdgeNum(); k++)
	{
		if (k != edgeIdx)
		{
			emx::real_data Rmk = sqrt(dot(Object->GetEdge(edgeIdx)->MidPoint - Object->GetEdge(k)->MidPoint,Object->GetEdge(edgeIdx)->MidPoint - Object->GetEdge(k)->MidPoint));
			emx::real_data t = j*Delta - Rmk/emx::C;
			if (t<0) continue;
			int floor = int(t/Delta);
			int ceiling = floor + 1;
			value += (T(floor, t) * I(k,floor) + T(ceiling, t) * I(k,ceiling)) * K(edgeIdx,k);
		}
	}

	return value;
}